### CREATED BY JESSICA SCHOENEHRR AND NICK WATERBURY ###
### Replication files for "Confessions at the Supreme Court" ###
### File 4 - Replication of MLM Models and Figures ###

### data sources:
### - Confessions of error originally collected by Bruhl (2009)
### 		- Updated by authors through the 2014 term
### - Politicization measure from Patrick Wohlfarth (2009)
###		- Updated by the authors through the 2014 term
### - sgJusticeDistance from the Judicial Common Space (Epstein, Martin, Segal, and Westerland)
###		- http://epstein.wustl.edu/research/JCS.html
### - CSI from Collins and Cooper
###		- https://dataverse.harvard.edu/dataset.xhtml?persistentId=doi:10.7910/DVN/UR2KYE
### - SG list of cases and SG as petitioner from Black and Owens (2012)
###		- updated by the authors per the OSG website
###		- https://www.justice.gov/osg/supreme-court-briefs
### - all other data from the Supreme Court DataBase (Spaeth, Epstein, et al.)
### 		- http://scdb.wustl.edu/index.php

library(arm)
library(compactr)
library(stargazer)
library(haven)

library(miceadds)
library(faraway)

library(tidyverse)
library(ggplot2)
library(ggthemes)
library(gridExtra)

##################
##################
##################
### MERITS MLM ###
##################
##################
##################

data1 <- read_dta("MeritsJusticeCentered20201026.dta")

data1$justice <- as.factor(data1$justice)
data1$votesWithSG <- as.factor(data1$votesWithSG)

### need to do this for PPs ###
data1$`92` <- ifelse(data1$justice == 92, 1, 0)
data1$`94` <- ifelse(data1$justice == 94, 1, 0)
data1$`95` <- ifelse(data1$justice == 95, 1, 0)
data1$`98` <- ifelse(data1$justice == 98, 1, 0)
data1$`99` <- ifelse(data1$justice == 99, 1, 0)
data1$`100` <- ifelse(data1$justice == 100, 1, 0)
data1$`101` <- ifelse(data1$justice == 101, 1, 0)
data1$`102` <- ifelse(data1$justice == 102, 1, 0)
data1$`103` <- ifelse(data1$justice == 103, 1, 0)
data1$`104` <- ifelse(data1$justice == 104, 1, 0)
data1$`105` <- ifelse(data1$justice == 105, 1, 0)
data1$`106` <- ifelse(data1$justice == 106, 1, 0)
data1$`107` <- ifelse(data1$justice == 107, 1, 0)
data1$`108` <- ifelse(data1$justice == 108, 1, 0)
data1$`109` <- ifelse(data1$justice == 109, 1, 0)
data1$`110` <- ifelse(data1$justice == 110, 1, 0)
data1$`111` <- ifelse(data1$justice == 111, 1, 0)
data1$`112` <- ifelse(data1$justice == 112, 1, 0)
data1$`113` <- ifelse(data1$justice == 113, 1, 0)
data1$`114` <- ifelse(data1$justice == 114, 1, 0)

justice <- data1[c(252:271)]

data1$`1` <- ifelse(data1$sgCounter == 1, 1, 0)
data1$`2` <- ifelse(data1$sgCounter == 2, 1, 0)
data1$`3` <- ifelse(data1$sgCounter == 3, 1, 0)
data1$`4` <- ifelse(data1$sgCounter == 4, 1, 0)
data1$`5` <- ifelse(data1$sgCounter == 5, 1, 0)
data1$`6` <- ifelse(data1$sgCounter == 6, 1, 0)
data1$`7` <- ifelse(data1$sgCounter == 7, 1, 0)
data1$`8` <- ifelse(data1$sgCounter == 8, 1, 0)
data1$`9` <- ifelse(data1$sgCounter == 9, 1, 0)
data1$`10` <- ifelse(data1$sgCounter == 10, 1, 0)
data1$`11` <- ifelse(data1$sgCounter == 11, 1, 0)
data1$`12` <- ifelse(data1$sgCounter == 12, 1, 0)
data1$`13` <- ifelse(data1$sgCounter == 13, 1, 0)
data1$`14` <- ifelse(data1$sgCounter == 14, 1, 0)
data1$`15` <- ifelse(data1$sgCounter == 15, 1, 0)
data1$`16` <- ifelse(data1$sgCounter == 16, 1, 0)
data1$`17` <- ifelse(data1$sgCounter == 17, 1, 0)
data1$`18` <- ifelse(data1$sgCounter == 18, 1, 0)
data1$`19` <- ifelse(data1$sgCounter == 19, 1, 0)

sgCounter <- data1[c(272:290)]

###########
### MLM ###
###########

model4MeritsFullMLM <- glmer(votesWithSG ~ policyChange
							+ noPolicyChange
							+ coeNoDetail
							+ coeCriminalCase
							+ coePastTally
							+ coeInPrevious
							+ politicization
							+ sgJusticeDist
							+ sgPetitioner
							+ constitutionalCase
							+ declareUncon
							+ CSI
							+ realOutlierSignal
							+ civilRightsCase
							+ (1 | justice)
							+ (1 | sgCounter),
							data = data1,
							family = binomial(link = logit))
							
summary(model4MeritsFullMLM)
logLik(model4MeritsFullMLM)
BIC(model4MeritsFullMLM)
nobs(model4MeritsFullMLM)

###
# policy change
###

set.seed(19870302)
n.draws <- 1000
sim.coef <- coef(sim(model4MeritsFullMLM, n.draws))
colnames(as.data.frame(sim.coef))
sim.coef <- as.data.frame(sim.coef)
colnames(sim.coef)

meritsPolicySeq <- seq(0, 1, length.out = 2)

pMeanMeritsPolicyNoCOE <- as.matrix(NA, nrow = length(meritsPolicySeq))
pLoMeritsPolicyNoCOE <- as.matrix(NA, nrow = length(meritsPolicySeq))
pHiMeritsPolicyNoCOE <- as.matrix(NA, nrow = length(meritsPolicySeq))

pMeanMeritsPolicyCOE <- as.matrix(NA, nrow = length(meritsPolicySeq))
pLoMeritsPolicyCOE <- as.matrix(NA, nrow = length(meritsPolicySeq))
pHiMeritsPolicyCOE <- as.matrix(NA, nrow = length(meritsPolicySeq))

pMeanMeritsPolicyCOEDiff <- as.matrix(NA, nrow = length(meritsPolicySeq))
pLoMeritsPolicyCOEDiff <- as.matrix(NA, nrow = length(meritsPolicySeq))
pHiMeritsPolicyCOEDiff <- as.matrix(NA, nrow = length(meritsPolicySeq))

for(i in 1:length(meritsPolicySeq)){
	predMeritsPolicyNoCOE <- cbind(1, 
							0,
							0,
							0,
							data1$coeCriminalCase,
							data1$coePastTally,
							data1$coeInPrevious,
							data1$politicization,
							data1$sgJusticeDist,
							data1$sgPetitioner,
							data1$constitutionalCase,
							data1$declareUncon,
							data1$CSI,
							data1$realOutlierSignal,
							data1$civilRightsCase,
							justice,
							sgCounter
							)
	predMeritsPolicyNoCOE <- as.matrix(predMeritsPolicyNoCOE)
	
	predMeritsPolicyCOE <- cbind(1, 
							1,
							0,
							0,
							data1$coeCriminalCase,
							data1$coePastTally,
							data1$coeInPrevious,
							data1$politicization,
							data1$sgJusticeDist,
							data1$sgPetitioner,
							data1$constitutionalCase,
							data1$declareUncon,
							data1$CSI,
							data1$realOutlierSignal,
							data1$civilRightsCase,
							justice,
							sgCounter
							)
	predMeritsPolicyCOE <- as.matrix(predMeritsPolicyCOE)
	
	ppMeritsPolicyNoCOE <- as.matrix(NA, row = n.draws)
	ppMeritsPolicyCOE <- as.matrix(NA, row = n.draws)
	ppMeritsPolicyCOEDiff <- as.matrix(NA, row = n.draws)
	
	for(j in 1:n.draws){
		xbMeritsPolicyNoCOE <- ilogit(predMeritsPolicyNoCOE %*% t(sim.coef[j, ]))
		ppMeritsPolicyNoCOE[j] <- mean(xbMeritsPolicyNoCOE)
		xbMeritsPolicyCOE <- ilogit(predMeritsPolicyCOE %*% t(sim.coef[j, ]))
		ppMeritsPolicyCOE[j] <- mean(xbMeritsPolicyCOE)
		ppMeritsPolicyCOEDiff[j] <- mean(xbMeritsPolicyCOE - xbMeritsPolicyNoCOE)
	}
	
	pMeanMeritsPolicyNoCOE <- mean(ppMeritsPolicyNoCOE)
	pLoMeritsPolicyNoCOE <- quantile(ppMeritsPolicyNoCOE, 0.025, na.rm = TRUE)
	pHiMeritsPolicyNoCOE <- quantile(ppMeritsPolicyNoCOE, 0.975, na.rm = TRUE)
	
	pMeanMeritsPolicyCOE <- mean(ppMeritsPolicyCOE)
	pLoMeritsPolicyCOE <- quantile(ppMeritsPolicyCOE, 0.025, na.rm = TRUE)
	pHiMeritsPolicyCOE <- quantile(ppMeritsPolicyCOE, 0.975, na.rm = TRUE)
	
	pMeanMeritsPolicyCOEDiff <- mean(ppMeritsPolicyCOEDiff)
	pLoMeritsPolicyCOEDiff <- quantile(ppMeritsPolicyCOEDiff, 0.025, na.rm = TRUE)
	pHiMeritsPolicyCOEDiff <- quantile(ppMeritsPolicyCOEDiff, 0.975, na.rm = TRUE)
}


pMeanMeritsPolicyNoCOE
pLoMeritsPolicyNoCOE
pHiMeritsPolicyNoCOE
nameMeritsPolicyNoCOE <- "No Confession of Error"

pMeanMeritsPolicyCOE
pLoMeritsPolicyCOE
pHiMeritsPolicyCOE
nameMeritsPolicyCOE <- "Policy Change"

pMeanMeritsPolicyCOEDiff
pLoMeritsPolicyCOEDiff
pHiMeritsPolicyCOEDiff

meritsBaseline <- data.frame(nameMeritsPolicyNoCOE, pLoMeritsPolicyNoCOE, pMeanMeritsPolicyNoCOE, pHiMeritsPolicyNoCOE)
colnames(meritsBaseline) <- c("div", "lb", "est", "ub")
rownames(meritsBaseline) <- c()

meritsPolicyCOE <- data.frame(nameMeritsPolicyCOE, pLoMeritsPolicyCOE, pMeanMeritsPolicyCOE, pHiMeritsPolicyCOE)
colnames(meritsPolicyCOE) <- c("div", "lb", "est", "ub")
rownames(meritsPolicyCOE) <- c()

###
# noPolicyChange
###

set.seed(19870302)
n.draws <- 1000
sim.coef <- coef(sim(model4MeritsFullMLM, n.draws))
colnames(as.data.frame(sim.coef))
sim.coef <- as.data.frame(sim.coef)
colnames(sim.coef)

meritsNoPolicySeq <- seq(0, 1, length.out = 2)

pMeanMeritsNoPolicyNoCOE <- as.matrix(NA, nrow = length(meritsNoPolicySeq))
pLoMeritsNoPolicyNoCOE <- as.matrix(NA, nrow = length(meritsNoPolicySeq))
pHiMeritsNoPolicyNoCOE <- as.matrix(NA, nrow = length(meritsNoPolicySeq))

pMeanMeritsNoPolicyCOE <- as.matrix(NA, nrow = length(meritsNoPolicySeq))
pLoMeritsNoPolicyCOE <- as.matrix(NA, nrow = length(meritsNoPolicySeq))
pHiMeritsNoPolicyCOE <- as.matrix(NA, nrow = length(meritsNoPolicySeq))

pMeanMeritsNoPolicyCOEDiff <- as.matrix(NA, nrow = length(meritsNoPolicySeq))
pLoMeritsNoPolicyCOEDiff <- as.matrix(NA, nrow = length(meritsNoPolicySeq))
pHiMeritsNoPolicyCOEDiff <- as.matrix(NA, nrow = length(meritsNoPolicySeq))

for(i in 1:length(meritsNoPolicySeq)){
	predMeritsNoPolicyNoCOE <- cbind(1, 
							0,
							0,
							0,
							data1$coeCriminalCase,
							data1$coePastTally,
							data1$coeInPrevious,
							data1$politicization,
							data1$sgJusticeDist,
							data1$sgPetitioner,
							data1$constitutionalCase,
							data1$declareUncon,
							data1$CSI,
							data1$realOutlierSignal,
							data1$civilRightsCase,
							justice,
							sgCounter
							)
	predMeritsNoPolicyNoCOE <- as.matrix(predMeritsNoPolicyNoCOE)
	
	predMeritsNoPolicyCOE <- cbind(1, 
							0,
							1,
							0,
							data1$coeCriminalCase,
							data1$coePastTally,
							data1$coeInPrevious,
							data1$politicization,
							data1$sgJusticeDist,
							data1$sgPetitioner,
							data1$constitutionalCase,
							data1$declareUncon,
							data1$CSI,
							data1$realOutlierSignal,
							data1$civilRightsCase,
							justice,
							sgCounter
							)
	predMeritsNoPolicyCOE <- as.matrix(predMeritsNoPolicyCOE)
	
	ppMeritsNoPolicyNoCOE <- as.matrix(NA, row = n.draws)
	ppMeritsNoPolicyCOE <- as.matrix(NA, row = n.draws)
	ppMeritsNoPolicyCOEDiff <- as.matrix(NA, row = n.draws)
	
	for(j in 1:n.draws){
		xbMeritsNoPolicyNoCOE <- ilogit(predMeritsNoPolicyNoCOE %*% t(sim.coef[j, ]))
		ppMeritsNoPolicyNoCOE[j] <- mean(xbMeritsNoPolicyNoCOE)
		xbMeritsNoPolicyCOE <- ilogit(predMeritsNoPolicyCOE %*% t(sim.coef[j, ]))
		ppMeritsNoPolicyCOE[j] <- mean(xbMeritsNoPolicyCOE)
		ppMeritsNoPolicyCOEDiff[j] <- mean(xbMeritsNoPolicyCOE - xbMeritsNoPolicyNoCOE)
	}
	
	pMeanMeritsNoPolicyNoCOE <- mean(ppMeritsNoPolicyNoCOE)
	pLoMeritsNoPolicyNoCOE <- quantile(ppMeritsNoPolicyNoCOE, 0.025, na.rm = TRUE)
	pHiMeritsNoPolicyNoCOE <- quantile(ppMeritsNoPolicyNoCOE, 0.975, na.rm = TRUE)
	
	pMeanMeritsNoPolicyCOE <- mean(ppMeritsNoPolicyCOE)
	pLoMeritsNoPolicyCOE <- quantile(ppMeritsNoPolicyCOE, 0.025, na.rm = TRUE)
	pHiMeritsNoPolicyCOE <- quantile(ppMeritsNoPolicyCOE, 0.975, na.rm = TRUE)
	
	pMeanMeritsNoPolicyCOEDiff <- mean(ppMeritsNoPolicyCOEDiff)
	pLoMeritsNoPolicyCOEDiff <- quantile(ppMeritsNoPolicyCOEDiff, 0.025, na.rm = TRUE)
	pHiMeritsNoPolicyCOEDiff <- quantile(ppMeritsNoPolicyCOEDiff, 0.975, na.rm = TRUE)
}


pMeanMeritsNoPolicyNoCOE
pLoMeritsNoPolicyNoCOE
pHiMeritsNoPolicyNoCOE
nameMeritsNoPolicyNoCOE <- "No Confession of Error"

pMeanMeritsNoPolicyCOE
pLoMeritsNoPolicyCOE
pHiMeritsNoPolicyCOE
nameMeritsNoPolicyCOE <- "Other Type of Confession"

pMeanMeritsNoPolicyCOEDiff
pLoMeritsNoPolicyCOEDiff
pHiMeritsNoPolicyCOEDiff

meritsNoPolicyCOE <- data.frame(nameMeritsNoPolicyCOE, pLoMeritsNoPolicyCOE, pMeanMeritsNoPolicyCOE, pHiMeritsNoPolicyCOE)
colnames(meritsNoPolicyCOE) <- c("div", "lb", "est", "ub")
rownames(meritsNoPolicyCOE) <- c()

meritsData <- rbind(meritsBaseline, meritsPolicyCOE, meritsNoPolicyCOE)

meritsData$div <- factor(meritsData$div, levels = c("No Confession of Error", "Policy Change", "Other Type of Confession"))

meritsFigure <- (ggplot(meritsData, aes(y = est, x = div, ymin = lb, ymax = ub)) +
	geom_pointrange(fatten = 1, size = 1) +
  	theme_light() +
 	ggtitle("SG as Party on the Merits") +
  	scale_x_discrete("", labels = c("No\nConfession\nof Error", "Policy\nChange", "Admission\nof\nMistake")) +
  	scale_y_continuous(limits = c(.50, .75), "Probability the Justice Votes with the OSG", breaks = seq(.50, .75, 0.02)) +
  	theme(aspect.ratio  = 1, plot.title = element_text(hjust = 0.5), panel.grid.minor.y = element_blank(), panel.grid.major.x = element_blank())
)

meritsFigure

##################
##################
##################
### VACB MODEL ###
##################
##################
##################

data3 <- read_dta("VoluntaryACBJusticeCentered20201026.dta")

data3$justice <- as.factor(data3$justice)
data3$votesWithSG <- as.factor(data3$votesWithSG)

### need to do this for PPs ###
data3$`92` <- ifelse(data3$justice == 92, 1, 0)
data3$`94` <- ifelse(data3$justice == 94, 1, 0)
data3$`95` <- ifelse(data3$justice == 95, 1, 0)
data3$`98` <- ifelse(data3$justice == 98, 1, 0)
data3$`99` <- ifelse(data3$justice == 99, 1, 0)
data3$`100` <- ifelse(data3$justice == 100, 1, 0)
data3$`101` <- ifelse(data3$justice == 101, 1, 0)
data3$`102` <- ifelse(data3$justice == 102, 1, 0)
data3$`103` <- ifelse(data3$justice == 103, 1, 0)
data3$`104` <- ifelse(data3$justice == 104, 1, 0)
data3$`105` <- ifelse(data3$justice == 105, 1, 0)
data3$`106` <- ifelse(data3$justice == 106, 1, 0)
data3$`107` <- ifelse(data3$justice == 107, 1, 0)
data3$`108` <- ifelse(data3$justice == 108, 1, 0)
data3$`109` <- ifelse(data3$justice == 109, 1, 0)
data3$`110` <- ifelse(data3$justice == 110, 1, 0)
data3$`111` <- ifelse(data3$justice == 111, 1, 0)
data3$`112` <- ifelse(data3$justice == 112, 1, 0)
data3$`113` <- ifelse(data3$justice == 113, 1, 0)
data3$`114` <- ifelse(data3$justice == 114, 1, 0)

justice3 <- data3[c(272:291)]

data3$`1` <- ifelse(data3$sgCounter == 1, 1, 0)
data3$`2` <- ifelse(data3$sgCounter == 2, 1, 0)
data3$`3` <- ifelse(data3$sgCounter == 3, 1, 0)
data3$`4` <- ifelse(data3$sgCounter == 4, 1, 0)
data3$`5` <- ifelse(data3$sgCounter == 5, 1, 0)
data3$`6` <- ifelse(data3$sgCounter == 6, 1, 0)
data3$`7` <- ifelse(data3$sgCounter == 7, 1, 0)
data3$`8` <- ifelse(data3$sgCounter == 8, 1, 0)
data3$`9` <- ifelse(data3$sgCounter == 9, 1, 0)
data3$`10` <- ifelse(data3$sgCounter == 10, 1, 0)
data3$`11` <- ifelse(data3$sgCounter == 11, 1, 0)
data3$`12` <- ifelse(data3$sgCounter == 12, 1, 0)
data3$`13` <- ifelse(data3$sgCounter == 13, 1, 0)
data3$`14` <- ifelse(data3$sgCounter == 14, 1, 0)
data3$`15` <- ifelse(data3$sgCounter == 15, 1, 0)
data3$`16` <- ifelse(data3$sgCounter == 16, 1, 0)
data3$`17` <- ifelse(data3$sgCounter == 17, 1, 0)
data3$`18` <- ifelse(data3$sgCounter == 18, 1, 0)
data3$`19` <- ifelse(data3$sgCounter == 19, 1, 0)

sgCounter3 <- data3[c(292:310)]


###########
### MLM ###
###########

model4VACBFullMLM <- glmer(votesWithSG ~ policyChange
							+ noPolicyChange
							+ coeNoDetail
							+ coeCriminalCase
							+ coePastTally
							+ coeInPrevious
							+ politicization
							+ sgJusticeDist
							+ sgPetitioner
							+ constitutionalCase
							+ declareUncon
							+ CSI
							+ realOutlierSignal
							+ civilRightsCase
							+ (1 | justice)
							+ (1 | sgCounter),
							data = data3,
							family = binomial(link = logit))
							
summary(model4VACBFullMLM)
logLik(model4VACBFullMLM)
BIC(model4VACBFullMLM)
nobs(model4VACBFullMLM)

###
# policy change
###

set.seed(19870302)
n.draws <- 1000
sim.coef <- coef(sim(model4VACBFullMLM, n.draws))
colnames(as.data.frame(sim.coef))
sim.coef <- as.data.frame(sim.coef)
colnames(sim.coef)

VACBPolicySeq <- seq(0, 1, length.out = 2)

pMeanVACBPolicyNoCOE <- as.matrix(NA, nrow = length(VACBPolicySeq))
pLoVACBPolicyNoCOE <- as.matrix(NA, nrow = length(VACBPolicySeq))
pHiVACBPolicyNoCOE <- as.matrix(NA, nrow = length(VACBPolicySeq))

pMeanVACBPolicyCOE <- as.matrix(NA, nrow = length(VACBPolicySeq))
pLoVACBPolicyCOE <- as.matrix(NA, nrow = length(VACBPolicySeq))
pHiVACBPolicyCOE <- as.matrix(NA, nrow = length(VACBPolicySeq))

pMeanMeritsPolicyCOEDiff <- as.matrix(NA, nrow = length(VACBPolicySeq))
pLoMeritsPolicyCOEDiff <- as.matrix(NA, nrow = length(VACBPolicySeq))
pHiMeritsPolicyCOEDiff <- as.matrix(NA, nrow = length(VACBPolicySeq))

for(i in 1:length(VACBPolicySeq)){
	predVACBPolicyNoCOE <- cbind(1, 
							0,
							0,
							0,
							data3$coeCriminalCase,
							data3$coePastTally,
							data3$coeInPrevious,
							data3$politicization,
							data3$sgJusticeDist,
							data3$sgPetitioner,
							data3$constitutionalCase,
							data3$declareUncon,
							data3$CSI,
							data3$realOutlierSignal,
							data3$civilRightsCase,
							#0 * data3$politicization,
							justice3,
							sgCounter3
							)
	predVACBPolicyNoCOE <- as.matrix(predVACBPolicyNoCOE)
	
	predVACBPolicyCOE <- cbind(1, 
							1,
							0,
							0,
							data3$coeCriminalCase,
							data3$coePastTally,
							data3$coeInPrevious,
							data3$politicization,
							data3$sgJusticeDist,
							data3$sgPetitioner,
							data3$constitutionalCase,
							data3$declareUncon,
							data3$CSI,
							data3$realOutlierSignal,
							data3$civilRightsCase,
							#1 * data3$politicization,
							justice3,
							sgCounter3
							)
	predVACBPolicyCOE <- as.matrix(predVACBPolicyCOE)
	
	ppVACBPolicyNoCOE <- as.matrix(NA, row = n.draws)
	ppVACBPolicyCOE <- as.matrix(NA, row = n.draws)
	ppVACBPolicyCOEDiff <- as.matrix(NA, row = n.draws)
	
	for(j in 1:n.draws){
		xbVACBPolicyNoCOE <- ilogit(predVACBPolicyNoCOE %*% t(sim.coef[j, ]))
		ppVACBPolicyNoCOE[j] <- mean(xbVACBPolicyNoCOE)
		xbVACBPolicyCOE <- ilogit(predVACBPolicyCOE %*% t(sim.coef[j, ]))
		ppVACBPolicyCOE[j] <- mean(xbVACBPolicyCOE)
		ppVACBPolicyCOEDiff[j] <- mean(xbVACBPolicyCOE - xbVACBPolicyNoCOE)
	}
	
	pMeanVACBPolicyNoCOE <- mean(ppVACBPolicyNoCOE)
	pLoVACBPolicyNoCOE <- quantile(ppVACBPolicyNoCOE, 0.025, na.rm = TRUE)
	pHiVACBPolicyNoCOE <- quantile(ppVACBPolicyNoCOE, 0.975, na.rm = TRUE)
	
	pMeanVACBPolicyCOE <- mean(ppVACBPolicyCOE)
	pLoVACBPolicyCOE <- quantile(ppVACBPolicyCOE, 0.025, na.rm = TRUE)
	pHiVACBPolicyCOE <- quantile(ppVACBPolicyCOE, 0.975, na.rm = TRUE)
	
	pMeanVACBPolicyCOEDiff <- mean(ppVACBPolicyCOEDiff)
	pLoVACBPolicyCOEDiff <- quantile(ppVACBPolicyCOEDiff, 0.025, na.rm = TRUE)
	pHiVACBPolicyCOEDiff <- quantile(ppVACBPolicyCOEDiff, 0.975, na.rm = TRUE)
}


pMeanVACBPolicyNoCOE
pLoVACBPolicyNoCOE
pHiVACBPolicyNoCOE
nameVACBPolicyNoCOE <- "No Confession of Error"

pMeanVACBPolicyCOE
pLoVACBPolicyCOE
pHiVACBPolicyCOE
nameVACBPolicyCOE <- "Policy Change"

pMeanVACBPolicyCOEDiff
pLoVACBPolicyCOEDiff
pHiVACBPolicyCOEDiff

VACBBaseline <- data.frame(nameVACBPolicyNoCOE, pLoVACBPolicyNoCOE, pMeanVACBPolicyNoCOE, pHiVACBPolicyNoCOE)
colnames(VACBBaseline) <- c("div", "lb", "est", "ub")
rownames(VACBBaseline) <- c()

VACBPolicyCOE <- data.frame(nameVACBPolicyCOE, pLoVACBPolicyCOE, pMeanVACBPolicyCOE, pHiVACBPolicyCOE)
colnames(VACBPolicyCOE) <- c("div", "lb", "est", "ub")
rownames(VACBPolicyCOE) <- c()

###
# noPolicyChange
###

set.seed(19870302)
n.draws <- 1000
sim.coef <- coef(sim(model4VACBFullMLM, n.draws))
colnames(as.data.frame(sim.coef))
sim.coef <- as.data.frame(sim.coef)
colnames(sim.coef)

VACBNoPolicySeq <- seq(0, 1, length.out = 2)

pMeanVACBNoPolicyNoCOE <- as.matrix(NA, nrow = length(VACBNoPolicySeq))
pLoVACBNoPolicyNoCOE <- as.matrix(NA, nrow = length(VACBNoPolicySeq))
pHiVACBNoPolicyNoCOE <- as.matrix(NA, nrow = length(VACBNoPolicySeq))

pMeanVACBNoPolicyCOE <- as.matrix(NA, nrow = length(VACBNoPolicySeq))
pLoVACBNoPolicyCOE <- as.matrix(NA, nrow = length(VACBNoPolicySeq))
pHiVACBNoPolicyCOE <- as.matrix(NA, nrow = length(VACBNoPolicySeq))

pMeanVACBNoPolicyCOEDiff <- as.matrix(NA, nrow = length(VACBNoPolicySeq))
pLoVACBNoPolicyCOEDiff <- as.matrix(NA, nrow = length(VACBNoPolicySeq))
pHiVACBNoPolicyCOEDiff <- as.matrix(NA, nrow = length(VACBNoPolicySeq))

for(i in 1:length(VACBNoPolicySeq)){
	predVACBNoPolicyNoCOE <- cbind(1, 
							0,
							0,
							0,
							data3$coeCriminalCase,
							data3$coePastTally,
							data3$coeInPrevious,
							data3$politicization,
							data3$sgJusticeDist,
							data3$sgPetitioner,
							data3$constitutionalCase,
							data3$declareUncon,
							data3$CSI,
							data3$realOutlierSignal,
							data3$civilRightsCase,
							#0 * data3$politicization,
							justice3,
							sgCounter3
							)
	predVACBNoPolicyNoCOE <- as.matrix(predVACBNoPolicyNoCOE)
	
	predVACBNoPolicyCOE <- cbind(1, 
							0,
							1,
							0,
							data3$coeCriminalCase,
							data3$coePastTally,
							data3$coeInPrevious,
							data3$politicization,
							data3$sgJusticeDist,
							data3$sgPetitioner,
							data3$constitutionalCase,
							data3$declareUncon,
							data3$CSI,
							data3$realOutlierSignal,
							data3$civilRightsCase,
							#0 * data3$politicization,
							justice3,
							sgCounter3
							)
	predVACBNoPolicyCOE <- as.matrix(predVACBNoPolicyCOE)
	
	ppVACBNoPolicyNoCOE <- as.matrix(NA, row = n.draws)
	ppVACBNoPolicyCOE <- as.matrix(NA, row = n.draws)
	ppVACBNoPolicyCOEDiff <- as.matrix(NA, row = n.draws)
	
	for(j in 1:n.draws){
		xbVACBNoPolicyNoCOE <- ilogit(predVACBNoPolicyNoCOE %*% t(sim.coef[j, ]))
		ppVACBNoPolicyNoCOE[j] <- mean(xbVACBNoPolicyNoCOE)
		xbVACBNoPolicyCOE <- ilogit(predVACBNoPolicyCOE %*% t(sim.coef[j, ]))
		ppVACBNoPolicyCOE[j] <- mean(xbVACBNoPolicyCOE)
		ppVACBNoPolicyCOEDiff[j] <- mean(xbVACBNoPolicyCOE - xbVACBNoPolicyNoCOE)
	}
	
	pMeanVACBNoPolicyNoCOE <- mean(ppVACBNoPolicyNoCOE)
	pLoVACBNoPolicyNoCOE <- quantile(ppVACBNoPolicyNoCOE, 0.025, na.rm = TRUE)
	pHiVACBNoPolicyNoCOE <- quantile(ppVACBNoPolicyNoCOE, 0.975, na.rm = TRUE)
	
	pMeanVACBNoPolicyCOE <- mean(ppVACBNoPolicyCOE)
	pLoVACBNoPolicyCOE <- quantile(ppVACBNoPolicyCOE, 0.025, na.rm = TRUE)
	pHiVACBNoPolicyCOE <- quantile(ppVACBNoPolicyCOE, 0.975, na.rm = TRUE)
	
	pMeanVACBNoPolicyCOEDiff <- mean(ppVACBNoPolicyCOEDiff)
	pLoVACBNoPolicyCOEDiff <- quantile(ppVACBNoPolicyCOEDiff, 0.025, na.rm = TRUE)
	pHiVACBNoPolicyCOEDiff <- quantile(ppVACBNoPolicyCOEDiff, 0.975, na.rm = TRUE)
}


pMeanVACBNoPolicyNoCOE
pLoVACBNoPolicyNoCOE
pHiVACBNoPolicyNoCOE
nameVACBNoPolicyNoCOE <- "No Confession of Error"

pMeanVACBNoPolicyCOE
pLoVACBNoPolicyCOE
pHiVACBNoPolicyCOE
nameVACBNoPolicyCOE <- "Other Type of Confession"

pMeanVACBNoPolicyCOEDiff
pLoVACBNoPolicyCOEDiff
pHiVACBNoPolicyCOEDiff

VACBNoPolicyCOE <- data.frame(nameVACBNoPolicyCOE, pLoVACBNoPolicyCOE, pMeanVACBNoPolicyCOE, pHiVACBNoPolicyCOE)
colnames(VACBNoPolicyCOE) <- c("div", "lb", "est", "ub")
rownames(VACBNoPolicyCOE) <- c()

VACBData <- rbind(VACBBaseline, VACBPolicyCOE, VACBNoPolicyCOE)

VACBData$div <- factor(VACBData$div, levels = c("No Confession of Error", "Policy Change", "Other Type of Confession"))

VACBFigure <- (ggplot(VACBData, aes(y = est, x = div, ymin = lb, ymax = ub)) +
	geom_pointrange(fatten = 1, size = 1) +
  	theme_light() +
 	ggtitle("SG as Voluntary Amicus") +
  	scale_x_discrete("", labels = c("No\nConfession\nof Error", "Policy\nChange", "Admission\nof\nMistake")) +
  	scale_y_continuous(limits = c(.50, .75), "Probability the Justice Sides with the OSG", breaks = seq(.50, .75, 0.02)) +
  	theme(aspect.ratio  = 1, plot.title = element_text(hjust = 0.5), panel.grid.minor.y = element_blank(), panel.grid.major.x = element_blank())
)

VACBFigure


### FIGURE ###

mlmFull <- grid.arrange(meritsFigure, VACBFigure, ncol = 2)




